The microscopic spectrum of the QCD Dirac operator with finite quark masses 
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We compute the microscopic spectrum of the QCD Dirac operator in the presence of dynamical 
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I. INTRODUCTION 



5h _ 

The low-lying eigenvalues of the QCD Dirac operator are of great importance for the understanding of the sponta- 
neous breaking of chiral symmetry in QCD. In Euclidean space, the Dirac operator reads ip — i$ + g(X a /2)^ a , where 
g is the coupling constant, A a are the generators of SU(iV)-color, and A° are the gauge fields. The corresponding 
£C) ' eigenvalue equation is ipipi — Xitpi with real eigenvalues A^. Because the Dirac operator anticommutes with 75, the 
eigenvalues either occur in pairs ±Aj or are zero. Let p(A) = (Ylj S(X — Aj))a be the eigenvalue density of the Dirac 
operator averaged over gauge field configurations. The Banks-Casher formula, (tpip) = irp(0)/V [0, relates the order 
parameter of chiral symmetry breaking, (V'V') 1 to the density of eigenvalues of the Dirac operator near zero virtuality. 
Here, it is essential that the thermodynamic limit be taken before the chiral limit. Note also that p(A) is normalized 
to the space-time volume V. In the real world chiral symmetry is spontaneously broken and (ipifi) 0- This implies 
that the low-lying Dirac eigenvalues must be spaced like ~ 1/V. Therefore, it is natural to magnify this region of the 
spectrum by a factor of V. This was first suggested by Shuryak and Verbaarschot § who introduced the so-called 
r~~- microscopic spectral density defined by 

o\ . 

2h ' Ps(z) = lim — p(— ) , (1) 

where E is the absolute value of the chiral condensate. Based on an analysis of sum rules derived by Leutwyler and 
Smilga H for the inverse powers of the Dirac eigenvalues in a finite volume, Shuryak and Verbaarschot conjectured 
that the microscopic spectral density should be a universal quantity which depends only on global symmetries, the 
number of flavors Nt, and the topological charge v. What is the meaning of universality in this statement? In 
principle, one could compute p s directly from full QCD. In practice, of course, the complexity of the QCD Lagrangian 
makes this impossible. However, if p s is universal, one can try to compute it in another theory which is simpler than 
QCD but has the same symmetries. The simplest theory in which such a calculation is possible is random-matrix 
theory (RMT). Whether the results obtained in RMT correctly describe the microscopic spectrum of full QCD is a 
question that can only be answered by empirical evidence, most importantly by large-scale lattice QCD simulations. 

Let us briefly recall how a chiral random-matrix model is constructed. In a chiral basis in Euclidean space, the 
Dirac operator has a block structure since it only couples states of opposite chirality. In RMT, the matrix of the 
Dirac operator is replaced by a random matrix W with suitable symmetry properties, 

(2) 

where we have also added a quark mass m. If W has N + v rows and N columns, then the Dirac matrix in (|^) 
(without the mass term) has N positive and N negative eigenvalues, respectively, and v zero modes (we assume v > 
and v <§; N). We can identify v with the topological charge and 2N + v w 2N with the volume V. In full QCD with 
Nf flavors, the weight function used in averaging over gauge field configurations contains the gluonic action in the 
form exp(— S s i) and Nf fermion determinants. In the random-matrix model, the gluonic part of the weight function 
is replaced by a convenient distribution of the random matrix W , usually a Gaussian distribution of the individual 
matrix elements (which are uncorrelated). Since the fermion determinants can be expressed as products over the Dirac 
eigenvalues they can also be taken into account in the random-matrix model. The symmetries of W are determined 
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by the anti-unitary symmetries of the Dirac operator. They were classified by Verbaarschot in Ref. 0). Depending on 
the number of colors and the representation of the fermions the matrix elements of W are real, complex, or quaternion 
real. The corresponding random-matrix ensembles are called chiral Gaussian orthogonal (chGOE), unitary (chGUE), 
and symplectic (chGSE) ensemble, respectively. The microscopic spectral density has been computed analytically in 
all three ensembles for all Nf and v in the chiral limit (^-Q]. 

Evidence in support of the conjecture that p s is a universal function has been accumulated in a number of recent 
studies which we list here. The moments of p s generate the Leutwyler-Smilga sum rules The result for p s 
is insensitive to the probability distribution of the random matrix g,^]. Lattice data for the valence quark mass 
dependence of the chiral condensate could be understood using the analytical expression for p s Jn],[ll| ■ The functional 
form of p s does not change at finite temperature fj^] . The analytical result for p s is found in the Hofstadter model 
for universal conductance fluctuations Jig ]. For an instanton liquid, p s shows quite good agreement with the random- 
matrix result . Recently, a high-statistics lattice calculation directly demonstrated agreement between the random- 
matrix results and lattice data, not only for p s but also for the distribution of the smallest eigenvalue and for the 
microscopic spectral two-point correlator (This study was done for gauge group SU(2) with staggered fermions 
in the quenched approximation.) We believe that the universality of p s is now firmly established. 

One of the most interesting and most important problems of present-day lattice simulations is the inclusion of 
dynamical fermions. This is essential for a realistic description of hadronic properties. Unfortunately, it is compu- 
tationally expensive to include the fermion determinants in a Monte-Carlo simulation. Moreover, the computational 
cost increases with lower quark mass, thus it is extremely difficult to take the chiral limit on the lattice. Any analytical 
insight into the distribution of the Dirac eigenvalues in the presence of massive dynamical quarks will therefore be 
helpful. It is the purpose of this work to extend the random-matrix results obtained in the chiral limit to the general 
case of massive quarks to be able to compare with actual lattice calculations. This issue was first discussed in Ref. [ pS[ . 
The most interesting case of QCD with 3 colors and fermions in the fundamental representation corresponds to the 
chGUE of RMT on which we concentrate here. In what range of the quark mass can we expect the random-matrix 
results to be in agreement with full QCD? Gasser and Leutwyler JlTj have shown that in the range 1/A <C L <C l/m n 
the dependence of the finite-volume QCD partition function on the quark mass m is determined solely by global 
symmetries (A is a typical hadronic scale, L is the linear size of the Euclidean box, and ~ V mA is the pion mass). 
The partition function was computed for equal quark masses using an effective Lagrangian in Ref. (^] and for the 
general case of different quark masses in the framework of chiral RMT in Ref. |L8| . The condition L <C 1 /m OT sets an 
upper limit of 1/VVA 2 on the quark mass. A lower limit of 1/(VA 3 ) on the quark mass is set by the requirement 
that the chiral condensate be non-zero in a finite volume. To obtain non-trivial results for the dependence of the 
various quantities we compute on m, the quark mass has to be rescaled by the same factor 1/(U£) as the eigenvalues, 
cf. Eq. (|lj). For large quark masses, our results reduce to those obtained in the quenched approximation. 

This paper is organized as follows. In Sec. [fl] we discuss the random-matrix model and construct the orthogonal 
polynomials and the partition function which are needed in the course of the calculation. Physical observables such as 



the microscopic spectral density and the distribution of the smallest eigenvalue are computed in Sec. III. Section [V 
is a short summary. 

After completion of this work we learned that very recently a similar calculation was performed simultaneously and 
independently by P.H. Damgaard and S. Nishigaki [|l9|. In this work, the authors also give a universality proof for 
more general distributions of the random matrix and construct spectral sum rules for the massive Dirac operator (see 
also Ref. Q ) but do not address the distribution of the smallest eigenvalue. Wherever the two manuscripts overlap 
the results are identical. 



II. RANDOM-MATRIX THEORY FOR FINITE QUARK MASSES 



In Sec. [I A , we formulate the chiral random-matrix model by intro ducing the probability density function and the 
partition function in the presence of finite quark masses. In S ec. II B, we construct the orthogonal polynomials of the 
model. An alternative construction is presented in Sec. II C which is based on a certain duality between models in 
different matrix spaces. 



A. Formulation of the model 



The statistical properties of the random-matrix model for the massive Dirac operator with Nf flavors and dynamical 
quark masses Mf (f = 1, . . . , Nf) which was introduced above are given by the probability density function 
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1 N f 

P { N S) {W,M) = — TT det^W^ + M 2 )exp(-AfE 2 trWW t ) (3) 

Z ( N Nf) (M) £ 

with M = (Mi, . . . , Mpf f ). Here, is a complex matrix of dimension AT. Apart from a Gaussian, this distribution 
contains Nf fermionic determinants. Eventually, we are mainly interested in the microscopic limit. Hence, it suffices 
to take into account a Gaussian in the total distribution (|^) since various authors have given proofs of universality 
in related problems. Note that even for nonzero topological charge v it is sufficient to consider only square matrices 
W in the calculation because of the duality between topology and flavor . A nonzero v is obtained by introducing 
v additional massless flavors at the end of the calculation. The normalization constant 

Z^ f \M) = / d[W] Yl det(WW f +Mf)exp(-NZ 2 trWW r ) (4) 
/=i 

depends on the quark masses and plays the role of a partition function. Since both of the above expressions depend on 
the matrices W only through invariant functions of WW^ , they can be reduced to integrals over radial coordinates. We 
write W = UAV with U G V(N), V € U(N)/J1 N (1), and A = diag(Ai, . . . , X N ), where the radial coordinates A, are 
restricted to the positive real axis. Thus, we have WW* = UA 2 U' f suggesting the introduction of new variables xi = A 2 . 
The measure transforms as d[W] = A 2 N (X)d[X]d^(U)d^{V), where X = diag(xi, . . . , xn), Ajv(X) = ni<j( x « — x j) 
is the Vandermonde determinant, and dfi is the Haar measure. After integrating over the unitary groups, we obtain 

1 N 

P { N Nf) (X,M) = A%(X) Y[wW( Xi ,M) (5) 

\M) i=l 

for the probability density and 

r N 

Z£'\M)= I d[X]A 2 N (X)l[wW(x i ,M) (6) 

i— 1 

for the partition function Zj^ f (M) which differs from (M) by the group volumes. The expression 

N f 

w {N f\x,M) = exp(-A^E 2 x) ]]_(x + Mf) (7) 

/=i 

will be referred to as the weight function. 

Because of a main result of RMT pi] , all spectral correlation functions R^ f \xi, . . . ,Xk,M) of our model can 

readily be written in terms of the polynomials p[^ f \x, M) which are orthogonal with respect to the weight function 

wW>(x,M), 

dxw( N f\x,M)p < h N f\x,M)p^\x,M) = S nm . (8) 

(n + l)4 Nf \M) 

The correlation functions are given as the determinant 

R ( k Nf) (xi,...,x k ,M) = det iK^ix^x^M)] (9) 

L J l,...,k 

with a kernel given by 

Kp\x,y,M) = yJw^f)(x,M)w^){y,M) ( " + p^\x, M)p^\y, M) . (10) 



With the help of the Christoffcl-Darboux formula the kernel can be expressed as 

z (N f ) 

(x, y, M) = N y/w&f)(x, M)w( N f)(y, M) 

Z N 

p^Hx, M)p$i\ (y, M) - P ™ (x, M)P$ f \v, M) 



x-y 

provided that the coefficient of the power x n in pi^ f \x, M) is unity. Eventually, we will take the microscopic limit 
as discussed in the Introduction 



(11) 
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B. Construction of orthogonal polynomials 



The polynomials pn* f \x, M) can be constructed by applying a standard formula |22j in the theory of orthogonal 
polynomials. In our case, it reads 



pi Nf) (x,M) = 



i r 

{Nfh - d[X}Al(X)Y[w( N fHx u M)( X - Xl ) 
Z n (M) J i=1 



(12) 



This is an integral over n variables Xi {i = 1, . . . , n). We notice that n replaces the dimension N of the matrices W 
in the partition function Z„ N (M) and in the Vandermonde determinant A n (X). Moreover, by construction, the 
coefficient in front of x n is unity. 

Fortunately, there is a further result, the Christoffcl formula [ p2[ , which allows us to immediately write down the 
polynomials explicitly: The weight function w^ Nf \x, M) is a product of the weight function w^°\x) and a polynomial 
whose zeros are the negative squared quark masses 



polynomials p^ f \x, M) can be written in terms of the polynomials which are orthogonal with respect to the weight 



Mj and which is nonnegative on the real axis. Therefore, the 



= exp(— NT, 2 x). These are precisely the Laguerre polynomials Ln' (NY, 2 x). The application of the Christoffcl 



(o), 



formula then yields 

pl Nf) (x,M) = 



1 1 (n + N f )\ 



l<£\ny?x) 

L ( °\-NE 2 M 2 ) 



L { n °UN^x) 
L^U-N^Ml 



L 



(0) 

»+ 

(0) 



Nf (NE 2 x) 



LZ Nf {~N^Ml 



(13) 



with a function Cn Nf \M) yet to be determined. In Sec. II C , we will present a very direct and highly convenient 



computation of the polynomials (13) including all normalization factors. Here we proceed as follows: Since the 



polynomial pif f \x, M) has to be invariant under permutations of the masses Mf and must have a finite limit when 
all quark masses are degenerate, a reasonable guess for the normalization function is 



C ( n Nf) (M) = det LWf^i-NYPM-j,) 



(0) 



fJ' = l,...,N f 



(14) 



see Sec. II C below, in particular Eq. (|24|). With this ansatz, the partition function Z„ / (M) is obtained 
by studying the case of Nf + 1 flavors. On the one hand, Eqs. (^) and ( |l2| ) yield p^ f (— Mjy +1 , M) = 
i)/Zn f \M), which, on the other hand, has to be equal to the right-hand side of Eq. (|l3|) 



-zL Nf+1) (M,M, 
for x = -Mf Jf+1 



N f + 1) 

Thus, we arrive at the recursion relation 



zL N ? +1 \M,M Nr 



l) Z { n Nl) {M) (n + N f )\ 



cL Nf+1) (M,M Nf 



C^\M) {N^Y^f UfUMj- M* f+1 ) 



(15) 



which implies 



n 



(n + f-l)\ C { n Ns) {M) 



L (iVS2)„+/-i A Nf {-M 2 ) 



(16) 



The correctness of these results can be verified using the orthogonality relation (g). In the microscopic limit, the 
partition function agrees with the result obtained in Ref. flS|| , see Sec. ID, Eq. (|27|). (Note that in Ref. |lj], the result 
is expressed in terms of derivatives of modified Bessel functions. Using the recurrence relations for Bessel functions 
|E3| and the properties of the determinant, one can show that the results are identical.) 
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C. A duality of matrix ensembles 



In this section, we present an alternative derivation of the orthogonal polynomials pn (x, M) and of the partition 
function Z^ 1 (M). It employs a certain kind of duality between matrix ensembles and is therefore of conceptual inter- 
est. On the practical side, this method advantageously yields all normalizations, particularly the function C^ ! \m) 
and the partition function Zn f \M), in a very direct and convenient way. We consider the integral 



r K 

X { *\a) = j d[W] exp(-atriyVK t ) det 

/=i 



— ICLf 



W 
-ia f 



(17) 



depending on K parameters aj (/ = 1, . 
partition function For N — n, K = 



, K). For a = NT?, K = Nf, and a f = Mf, the integral ([17]) is just the 



N f +!,«/= Mf (/ = !,..., N f ), and a? Nf 



+1 



up to the group volumes, equal to Zn' <f> (M)p ( ^ f> (x,M) as denned in Eq. fll2]). Thus, the function X)^ > (a) yields 



-x, the integral ( |17| ) is, 



all desired quantities. We define K 2A^-component vectors £/ (/ = 1, . 
(i = 1, . . . , 27V) as well as the 2iV.fr-component combination £ = (Ci, ■ 
Eq. (|l7|) are represented as the integral 



. , K) with complex Grassmannian entries Qi 
, Ck) T ■ The determinants in the integrand of 



K 



JJdet 



— Idf 



w 

-ia f 



(2tt) 2NK J d[Q exp ^ (l K 







IT' 




(18) 



with a = diag(ai, . . . ,a/{). In this form, the Gaussian average over the matrices W can easily be performed. We 
arrive at a four-fermion interaction which is then decoupled by a Hubbard-Stratonovitch transformation. In Ref. p4| , 
a very similar route was taken. In this study, however, a problem involving fermionic and bosonic determinants was 
solved which led to the introduction of supervectors of commuting and anticommuting variables. In the present case, 
the vectors Q have only anticommuting entries. Fortunately, this difference has little influence on the structure of 
the results. Instead of the supermatrices which were used in Ref. 24 for the Hubbard-Stratonovitch transformation, 
we have to employ ordinary matrices here. We arrive at 



X 



N 



(a) = ^— ^ / d[a] exp (— atr (a — a) (a^ — a)) det 



N f 

era 1 



(19) 



whe re a is an ordinary complex KxK matrix without further symmetries. Thus, we have mapped the original model of 
Sec. [I A in the form ([l7]) in the space of the NxN matrices W onto a dual model in the space of the KxK matrices a. 
This is highly convenient because, eventually, we wish to ta ke the limit N — > oo. Remarkably, these dual models differ 
with respect to rotation invariance. The model of Sec. II A is invariant under rotations in the sense that the integrand 
of the partition functions in Eqs. (Q) and (S) depends only on the radial coordinates of the matrices W, i.e., on the 
variables xi. This is not true for the model (Jl9|). The presence of the matrices a breaks this kind of rotation invariance. 
In Refs. [^D and (24), a procedure was developed to deal with models precisely of this kind. We introduce radial 
coordinates by a = usv with u e U(K), v e V(K)/V K (1), and s = diag(si, . . . , s K ). A gain, the radial coordinates 
Sf are restricted to the positive real axis. The measure transforms as d[a] = A 2 K (s 2 ) Il/Li Sfd[s]dfj,(u)diu,(v). The 
integral over the unitary groups is now nontrivial. The solution can be found in Ref. |2q|, 



^— ^ / dfi(u) / d/i(v) exp (—atr (usv — a)(u^stt^ — a)) 



1 det[7(s / ,a / >)] /J/=li 



A" 



.71"/ j j Kl A K (s 2 )A K (a 2 ) 

The entries of the determinant are given by 

■y(sf,aff) — 2a exp (— a(s 2 f + a 2 ,)) I (2asfa,f>) , 
where Iq is the modified Bessel function of the first kind and zeroth order. Collecting everything, we obtain 



K\ A K (a 2 ) 



K poo 

[] / dsfs} N+1 A K ( S 2 )det[-y(s f , 
f=i Jo 



*f')\f, f >.. 



.K 



(20) 



(21) 



(22) 



This integral can be solved in a straightforward way. We write Ak(s 2 ) 
determinant row by row. Using the representation 



det Is 



2(/'-l) 
/ 



]/,/'=!, ...,K and integrate this 
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2a / dw w 
Jo 



2n+2m+l 



exp(-a(z 2 +W 2 ))/ (2a W z) = L^-az*) 



for the Laguerre polynomials, we find 
X f N K \a) 



1 n'f^det 



A K (a 2 ) -11 a N+f 
v ; /=o 



-(o) 

J N+f- 



-1 ) 



/,/'=l,...,K 



This result immediately yields the orthogonal polynomials p < ^' f \x, M) and the partition function Z^ ! \m). 



III. CALCULATION OF OBSERVABLES 



(23) 



(24) 



T he mi croscopic spectral correlators, in particular the microscopic spectral density, are computed in Sec. Ill A. In 
Sec. [II B, we derive the distribution of the smallest eigenvalue. 



A. Microscopic spectral correlators and density 



To calculate the spectral correlations, we insert the polynomials pl^ !f \x 1 M) in Eq. (|ll| ) for the kernel 
Kj^ f (x,y, M). The microscopic limit is obtained according to Eq. ([!]) by rescaling energies and quark masses 
by 2NH. We define A = 2NHy/x and to/ — 2NHMf. In the limit N — » oo, we want both A and to/ to be of order 
unity. The quenched approximation then corresponds to taking to/ — > oo (or, of course, setting Nf = 0). In the limit 
to/ — > 0, our results must reduce to those obtained earlier in the chiral limit 0|. After rescaling the parameters one 

uses the asymptotic relation limjv^oo \z/N)/N a = z a/2 J a (2^/z), where J a denotes the Bessel function of the 
first kind and order a. To obtain the correct large- N limit one has to use the recurrence relations for the Laguerre 
polynomials (23) in Eq. (|l3|). The properties of the determinant can be used to simplify the expression. After summing 
up one uses 



N 



N ■ 



L 



N 



N ■ 



L 



N 



(a+l) 
N 



N 



N ■ 



a+l 



-(a) 
J N 



N 



(25) 



Note that the first, second, and third term on the right-hand side is of 0(1), 0(1/N), and 0(1/A 2 ), respectively. 
Using the properties of the determinant and collecting the terms of 0(1) carefully, one obtains the microscopic limit 
of (fn]). We arrive at 



K {N i ] {\ A', to) = lim 



2VAA' 



K 



(N f ) 



((A/2AE) 2 , (A'/2A£) 2 ,to/2AI]) 



N-T™ (2 AS) 2 N 
AA 7 B( N f\X, to) J2% Bf'\\', to) - flW(A', to) E£o Bf { \\, > 



A' 2 - A 2 



[C(^/)(to)] 2 Yf f U y / (A 2 + TO 2 )(A' 2 +m 2 ) 



(26) 



where the normalization C^- N f\m) is the microscopic limit of the normalization Cjy (M) of Eq. (fbjl) 



CW(m) = lim C^ V/) (TO/2A f S) 



det 



f,f>=l,...,Nf 



(27) 



The determinant B^ N f^ (A, to) in Eq. ( p6j ) is given by 



^o(A) Jo(-TOi) 

AJl(A) TOl/l( — TOl) 



Io(-m Nf ) 
mN f h{-m Nf ) 



X N fJ Nf (X) m^ f I Nf (~mi) ■■■ m N f f I Nf (-m Nf ) 



Nf; 



(28) 



G 



The determinant 
is constructed from B( N f'(X,m) by replacing the first column by (AJi(A), 



D 



(N f ) 



(A,m) 



replacing the (/+l)-th column by (m/Ji(— m/), . . . ,m,' I]s f +i{—mf)) T if / = 1 



relations can be obtained from the kernel ( |26| ) in analogy with Eq. (|9j 
density is obtained by taking the limit A' = A = z in (Eq). This yields 



,A Ar /+ 1 J JV/+1 (A)) T if / = and by 
, Nf . All fc-point spectral cor- 
Most importantly, the microscopic spectral 



pi j; (z,to) = - 



Nf jj 
f=0 U f 



(Nf) 



(z, to) - DW(z, to) E/£ B f" ( z > m ) 



7(Nf), 



[c( N /)(TO)] 2 n7ii( z2 + TO /) 



(29) 



where the determinant D^ N f\z,m) is constructed from B^ N ^ (z,m) by replacing the first column in ( p8| ) by 
(z -1 J_i(z), . . . , z^" 1 Jat / _ 1 (z)) t , Dj Nf \z,m) is constructed from D^ N ^(z,m) by replacing the (/+l)-th column 

by (TO / / 1 (-TO / ),...,mJ r/+1 / A r /+1 (-TO / )) T if /= l,...,N f , and D^ Nf \z, m) = BW) (z, to). 

We have verified that our results ( p(f ) and (^9|) reproduce the known results for the chiral limit and the quenched 
approximation if to/ — ► and to/ — > oo, respectively. In Refs. [ p4|j2^1 , the correlations of the Dirac operator at finite 
temperature were studied using the supersymmetry method. These investigations suggest that the expression ( p^ ) 
for the kernel can be simplified further and expressed in form of a single determinant. Indeed, such a more compact 
form was very recently obtained in Rcf. |19| . We write the result in the form 



K {N f\\,\',m) = 



'XX' 



1 



A' 2 - A 2 



(to) y / (A 2 + TO 2 )(A' 2 +to 2 ) 



Jo (A) 
AJi(A) 



Jo(A') 
A'Ji(A') 



^o(-TOl) 
TOiJi(-TOi) 



A^+V^+iCA) A'^ +1 J^ +1 (A') mf' +1 Jjv /+ i(-mi) 
The microscopic spectral density then reads 



Io(-m Nf ) 
m Nf Ii(-m Nf ) 



Nf+l T I \ 

jv> 1 N f +i{-m Nf ) 



(30) 



p s (z,to) = -~ 



2c( jv /)(TO)n74(^ 2 +^/) 

Jq(z) zJi(z) 



7 (-mi) 

TOlil(-TOl) 



z N fJ Nf (z) z N f +1 J Nf+1 (z) mf f+1 I Nf+ i(-mi) 



I (-m Nf ) 
mN f Ii(-m Nf ) 



m 



Nf 



1 lN f +i(—i r nN f ) 



(31) 



We have verified by direct calculation that the expressions (p6[) , ( j 30| ) and (29), (|3l]) are identical, respectively. A 
number of interesting special cases can be derived from (^9|) or (|3l|). (For Nf = 1 and v = 0, p s was first computed 
in Ref. However, like the authors of 19 1, we fail to agree with the result given in Eq. (80) of Rcf. Jl6| which does 
not reproduce the quenched result in the limit to — > oo.) As an example, we give the result for p s for two degenerate 
flavors of mass to and v = 0, 



„(2) 



(z,to) = -(J 2 (z) + J 2 (z)) 



2z 



[zJi(z)I Q (m) + mJ (z)h(m)f 



(z 2 + TO 2 ) 2 



Ii(m)-I((m) 



(32) 



This expression is plotted in Fig. [I] compared with the results obtained in the chiral limit and in the quenched 
approximation. Clearly, the presence of light flavors leads to a suppression of small eigenvalues as expected. We 
emphasize again that while we have worked with v = 0, the general case v > can easily be obtained by introducing 
v additional massless flavors. For completeness, we state the result for the microscopic spectral density for general v 
which is obtained by expanding the determinants in (|3l|), 

(Nf,v)( n 

p s (z,to) 



2c(Nf,"){rn) z ^ ri/=i^ 2 + ™/) 
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Z V ~ X J v -x{z) 

Z v J v {z) 



z v J v {z) 
z v+1 J v+1 {z) 



m\ +1 I v+ i{-mi) 



™>N f 1 Iu+i(-'m Nf ) 



z N f+ u J Nf+v (z) z N f + " +1 J Nj+lJ+ i(z) mf t+l/+1 I Nf+u+ i(-m 1 ) ... m^ f +u+1 I Nf+u+ i(-m Nf ) 



(33) 



with 



& N f")(m) = det \m v f t S ~ X I v+ f-i(-m r ) 



f,f>=l,...,N, 



(34) 



Note that the chiral limit follows immediately from ( |33| ) by the replacements Nf — and v — * Nf + v. We then 
obtain 



pi Nf ' v) (z,0) = | [J Nf+u {z) - J N/+u+1 (z)J Nf+u -i(zyj 



(35) 



as required ||. 



B. Distribution of the smallest eigenvalue 

In lattice simulations, the performance of an algorithm is frequently determined by the magnitude of the smallest 
eigenvalue. It is therefore of interest to compute an analytical result for the distribution of the smallest eigenvalue, 
p( N f> (A m i n , m), in the framework of RMT. This is the subject of this section. The RMT result should be universal and 
identical with that of full QCD under the same conditions as the microscopic spectral density. The only result that 
has been computed so far in the literature applies to the quenched approximation where p(°\\) — Aexp(— A 2 /4)/2 

©■ 

A well-known quantity in RMT is the so-called "hole probability" E(s\, S2) which is the probability that the interval 
(si, S2) is free of eigenvalues. For the smallest eigenvalue, we again move to the microscopic scale by rescaling energies 
and masses by 2iVE. For the interval (0, s) at the spectrum edge, we then have 



E(0,s) = / dxi ■ ■ ■ dxNP^ f \xi, . . . ,XN',rn) 



where the probability density function on the microscopic scale (with Xi = zf) is given by 



N 



Nj 



^ f) (x u . ..,x N ;m) = c { p\m)A%(x) [J e~^ iN J] (z, + mj) 



(36) 



(37) 



/=i 



cf. Eqs. (§) and (0). The factor c^ f \m) ensures that E(0,0) = 1 and should not be confused with the normalization 
constant defined in Eqs. ( |l4| ) and (|271). On the other hand, we have 

(38) 



E(0,s) = J dx mia P^{ 
Hence, the distribution of the smallest eigenvalue (x m i n = A^ in ) is given by 

P^(Amin,m) = 2A min F( A '^(.T min ,m) = -2A 
Performing the derivatives, we obtain 

- E'(0, s) = c^ f) (m)Ne^ 4N jj ( s + m}) 



dE{0, s) 



ds 



(39) 



/=i 



N-l 



/ d Xl ■ ■ ■ dxN-x^ N ^{x) J] e- x ^ iN { Xl - sf TJ ( Xi + m}) 

J 8 • 1 



Nf 

n 

/=i 



(40) 
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We now shift x% — > x% + s and finally obtain from 



Nf 

(A, m) = 2c<£ f) (m)N\e- x2 ^ f[ (A 2 + to 2 ) 

/=i 

,00 N-l N f 

x dx 1 ---dx N . 1 A 2 N _ 1 (x)Y[e- x ' /iN x 2 l Y[(x l + X 2 +m 2 f ), (41) 
Jo i=l /=i 

where we have written A instead of A m in for brevity. The integrals in Eq. ( |4l"| ) can be done by comparison with the 
case of the microscopic spectral density. In terms of the probability density function on the microscopic scale, p s is 
given by 



p s f '(z, to) = 2z ■ N / dxi ■ ■ ■ dxN-i]r N (xi, ■ ■ ■ ,xn-i, z ; to) 
Jo 

N f 

JNf), 



2c { „ f) (m)NzJ[(z 2 + m 2 ) 
f=l 

N-l N f 

(42) 

f=i 



,00 N-l Nf 

/ dx 1 ---dx N - 1 A 2 N _ 1 (x)l[e-*^ m (x l -z 2 ) 2 l[(x l +m 2 ) 
Jo i=i f=i 



in the large- N limit. Apart from prefactors in front of the integrals, we note that P( N ^(\,m) can be obtained by 
setting z = in the integrand of Eq. ( ]42] ) and by replacing to/ by (A 2 + to 2 ) 1 / 2 in the result for p s (but not in the 
normalization factors). Setting z = in the numerator determinant of Eq. ( |3l| ) makes the first two columns trivial, 
and the dimension of the matrix can be reduced to Nf by expanding the determinant. Including the prefactors and 
the proper normalization, we obtain 



(A 2 + to 2 ,) (/+1)/2 I f+1 (-(A 2 + to 2 ,) 1 / 2 ) 



det 

= r~ x2/i — ~. ' AfJ ' =1 --' Nf m->- 



C( A '/)(TO)n/ii(A 2 + mj) 



with C*W(m) given in Eq. (p7|). For Nf = 0, Eq. (|4^) reduces to the quenched result trivially. We also reproduce 
the quenched result by taking the limit rrif — * 00. The simplest nontrivial case is Nf = l,v = where 

pW(A,m) = ^ ^ + g) . (44) 
2 J (m) 

For two degenerate flavors of mass m and ^ = 0, we obtain 

In Fig. ||, we plot the distribution of the smallest eigenvalue for the same parameters as in Fig. |l|. We observe the 
same suppression of small eigenvalues in the presence of light dynamical quarks. Other cases can readily be derived 
from ( p3| ) . The most general case where v > can again be obtained by introducing v additional massless flavors and 
expanding the determinants in (|43|). As an example, we state the result for one massive flavor and v — 1, 



P (iv / =i,-i) (A)m) = ^ e -AV4_L_ [t j a (A)J 3 (t) - A / 2 (*)/3(A)] t=VIW . (46) 



IV. DISCUSSION 



We have computed the microscopic spectrum of the massive QCD Dirac operator, including the microscopic spectral 
density and the distribution of the smallest eigenvalue, for an arbitrary number of flavors, arbitrary quark masses, and 
arbitrary topological charge. Our results generalize the previously known results for the microscopic correlations in the 
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chiral limit j5j . The distribution of the smallest eigenvalue was previously known only in the quenched approximation. 
Wherever there is an overlap, our results agree with those obtained in the recent preprint by Damgaard and Nishigaki 
p9| . It would be very interesting to obtain the result (|3^) using the supersymmetry method as in Refs. |24j,^6|. Work 
in this direction is in progress. 

The quantities we have computed in this paper are conjectured to be universal, i.e., identical to those of QCD, 
for quark masses of the order of l/V. This conjecture can be verified straightforwardly in lattice calculations. One 
computes the eigenvalues of the Dirac matrix and constructs the microscopic spectral density and the distribution of 
the smallest eigenvalue by averaging over a sufficient number of independent gauge field configurations. The lattice 
data can then be compared to the analytical predictions of RMT. For pure SU(2) gauge theory without dynamical 
quarks, this was already done in Ref. |[15| , and the agreement was excellent. 

Our present calculation applies to full QCD with 3 colors in the continuum limit. In the context of lattice QCD, it 
applies to SU(3) with dynamical staggered fermions. The necessary lattice algorithms are available, e.g., the hybrid 
Monte-Carlo algorithm for generating unqucnchcd configurations and the Lanczos algorithm for computing eigenvalues 
of large sparse matrices. Therefore, we are confident that the results obtained in this paper will be verified in the 
very near future. 

Given that the analytical results obtained in this paper really describe QCD, one may ask what their relevance is 
in practice. We list here a few possible applications, (i) The parameter which appears in the microscopic quantities 
contains the absolute value of the chiral condensate, E. As we have already seen in quenched SU(2) |2§], the 
analytical information available from RMT can be used to determine £ and to extract its thermodynamic limit from 
very small lattices. Since otherwise one would have to perform an extrapolation, this seems to offer an interesting 
technological advantage, (ii) Lattice simulations with light dynamical quarks are very expensive, and extrapolations 
to the chiral limit are difficult. In RMT, we have computed exact analytical expressions which apply for very small 
quark masses, and the chiral limit can be taken analytically. One would, therefore, hope that lattice practitioners 
can use the analytical information in their extrapolations to the chiral limit for those quantities that depend on the 
small eigenvalues, (iii) A very interesting issue is topology. As one approaches the continuum limit on the lattice, 
one expects zero modes to appear which give rise to non-zero topological charge v. The microscopic quantities we 
consider are sensitive to topology, and we have computed them for arbitrary v. Hence, one can analyze lattice results 
in different sectors of topological charge. Whether or not topology is seen on present-day lattices is still a somewhat 
controversial issue, and for quenched SU(2) the lattice results of Ref. |l5| were consistent with v = 0. We hope that 
our results will prove useful in future analyses of topology, (iv) The performance of many lattice algorithms depends 
on the magnitude of the smallest eigenvalue. While RMT cannot predict this magnitude, it predicts the distribution of 
Amin over the ensemble of gauge field configurations. It would be interesting to see if experts in algorithm development 
can make use of the analytical information available from RMT. 
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FIG. 1. Microscopic spectral density for two degenerate flavors of rescaled mass m = 5 (solid line), in the chiral limit (short 
dashes), and in the quenched approximation (long dashes). 
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FIG. 2. Distribution of the smallest eigenvalue for two degenerate flavors of rescaled mass m = 5 (solid line), in the chiral 
limit (short dashes), and in the quenched approximation (long dashes). 
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